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ABSTRACT 

A computer code for solving the Reynolds averaged full Navier-Stokes equations has 
been developed and applied using sheared H-type grids. The Baldwin-Lomax eddy-viscosity 
model is used for turbulence closure. The integration in time is based on an explicit 
four-stage Runge-Kutta scheme. Local time stepping, variable coefficient implicit residual 
smoothing, and a full multigrid method have been implemented to accelerate steady state 
calculations. Comparisons with experimental data show that the code is an accurate vis- 
cous solver and can give very good blade-to-blade predictions for engineering applications 
in less than 100 multigrid cycles on the finest mesh. 


*This research was supported by the National Aeronautics and Space Administration under Contract 
No. NAS1-18107 while the first author was in residence at the Institute for Computer Applications in Science 
and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 







INTRODUCTION 

The capability to correctly predict the various flow conditions that can occur in tur- 
bomachinery blade passages is a basic step for the design and the improvement of the 
modern compressors and turbines. As indicated by the research of the past few years, 
flow solutions can be obtained rapidly by solving the potential equation, but this equation 
does not allow for vorticity and entropy changes related to shock waves. Recently, many 
numerical methods have been developed which solve the Euler equations. However, an 
inviscid analysis requires some assumption about the circulation to get an unique solution, 
and the problem may become ambiguous when dealing with a rounded trailing edge. In 
addition, heat transfer and boundary-layer thickening, that can lead to separation and hot 
spots, cannot be represented in an inviscid solution. 

By coupling a boundary-layer model with an inviscid procedure one can remove some 
of the previous limitations, but such an approach is difficult to extend to three-dimensional 
flows. The improvement in the numerical schemes and in the performance of modern com- 
puters suggests investigating the solution of the Reynolds averaged Navier-Stokes equa- 
tions. Although significant progress has been made in the development of Navier-Stokes 
solvers for cascade flows, as evident from [1-3], much work remains before computer codes 
capable of efficiently solving the wide range of blade geometries and flow conditions that 
can be encountered in the design process are available. 

The objective of the present work is to develop a fast and accurate method for the 
computation of viscous cascade flows. A four stage Runge-Kutta scheme with acceleration 
techniques for steady-state solutions is used to advance in time the Reynolds averaged 
Navier-Stokes equations. Numerical efficiency has been achieved through local time step- 
ping, variable coefficient implicit residual smoothing, and a full multigrid method. Sheared 
H-type grids are used to discretize the computational domain. These grids are easy to 
generate algebraically and to extend to three-dimensional problems, and they generally 
provide sufficient resolution for engineering purposes. 

In this paper, numerical results are presented for a bicircular arc cascade, for a typical 
gas turbine high loaded cascade with thick, rounded trailing edge, and for a bump in a 
channel. The bicircular profile is studied for the case of subsonic, laminar, separated flow 
and for the case of transonic, turbulent flow with separation at the foot of the shock and 
near the trailing edge. For the gas turbine blade, subsonic and transonic conditions are 
considered, and solutions are compared to experimental data. Finally, the capability of the 
present prediction technique to compute supersonic flows is demonstrated with the bump 
in a channel problem. For all these cases, convergence histories are shown. 
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GOVERNING EQUATIONS 

Let p,u,v,p,T,E, and H denote respectively the density, velocity components in the x 
and y cartesian directions, pressure, temperature, specific total internal energy, and specific 
total enthalpy. The two-dimensional, unsteady Navier-Stokes equations, neglecting body 


forces and heat sources, can be written in conservative form in a cartesian coordinate 
system as, 

dU dF dG dR dS 

dt dx dy dx dy V ’ 

where 


and 



p 


pu 


pv 

u = 

pu 

; F = 

pu 2 + p 

; G = 

puv 


pv 


puv 


pv 2 + p 


pE 


puH 


pvH 



0 " 


0 

R = 

T X z 

s = 

Tyx 


T X y 


Tyy 


r 4 


. S * 


r X x = (A + 2p)u x + Av„ 
T w = (A + + A u x 

T X y = T yX p(u y + Vx) 

r 4 = UT XX + VT XV + kT x 

S 4 = UTy X + VTyy + kTy. 


( 2 ) 

( 3 ) 
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The pressure is obtained from the equation of state, 


p = pRT. 


( 5 ) 


According to the Stokes hypothesis A is taken to be —2/3p and a power law is used to 
determine the molecular coefficient of viscosity p. The effect of turbulence is taken into 
account by using the eddy-viscosity hypothesis. That is, the molecular viscosity p and the 
molecular thermal conductivity k are replaced with, 



where c p is the specific heat at constant pressure, Pr is the Prandtl number, and the 
subscripts £ and t refer to laminar and turbulent. The turbulence quantities nt and Pr t 
are computed using the two-layer algebraic model of Baldwin and Lomax [4]. 

The Euler equations can be easily obtained by neglecting the right hand side of equation 

( 1 ). 


SPATIAL DISCRETIZATION 

A finite-volume approach is applied to discretize the equations of motion. The com- 
putational domain is divided into quadrilateral cells, fixed in time, and for each cell the 
governing equations are written in integral form as follows, 



Vixiv + jjFdy - Git) = fjRdy - Sdx) 


( 8 ) 


where fi is a generic cell and dfl its boundary. 

A cell-centered discretization is used, and the line integral of (8) is approximated with 
the midpoint rule. The convective fluxes at the cell faces are obtained by simple averaging 
of adjacent cell-center values of the dependent variables (see [5]). The spatial derivatives 
necessary to compute the viscous terms are evaluated by means of Green’s theorem. For 
example, consider the arbitrary cell fi (EFGH) in Fig. 1. The contributions u x and u y to 
the viscous flux across each cell face (i.e., GH) are determined using 



( 9 ) 


where fi' is an auxiliary cell (CGNH for face GH in Fig. l). In a similar manner, v x and v v 
are computed. Additional details for finite-volume treatment of viscous stresses and heat 
conduction terms are found in [5] and [6]. 


BOUNDARY CONDITIONS 

It is well-known that the reflective behavior of numerical boundary conditions can sig- 
nificantly influence the convergence rate of a marching scheme. Thus, boundary conditions 
must be chosen carefully. According to the theory of characteristics, three quantities are 
specified at a subsonic inflow boundary. In the present work, the total enthalpy, the total 
pressure, and flow angle are specified if the axial flow is subsonic; the outcoming Riemann 
invariant R~ is taken from the interior. For supersonic axial inlet all the flow variables 
are computed from the prescribed inlet conditions. At a subsonic axial outlet, the static 
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pressure corresponding to the desired isentropic Mach number is prescribed, while the out- 
going Riemann invariant R + , the total enthalpy, and the component of the velocity parallel 
to the boundary are extrapolated from the interior. If the axial outlet is supersonic all 
the variables are extrapolated. On the solid wall, the only contribution to the convective 
fluxes is the pressure, which is obtained by extrapolation from the interior of the flow field. 
For viscous flow calculations, the velocity is set to zero on the solid surface, which is taken 
to be adiabatic in the present work. 


ARTIFICIAL DISSIPATION 

For inviscid flows the finite-volume scheme does not contain any dissipative terms. In 
order to prevent odd-even point decoupling and oscillations near shock waves or stagnation 
points, artificial dissipation terms are added to the governing discrete equations. For 
viscous flows dissipative properties are present due to diffusive terms; however, away from 
shear layer regions (i.e., boundary layers, wakes) the physical dissipation may not be 
sufficient to guarantee stability, especially in the case of stretched meshes. Although shock 
waves can be captured by defining the shock structure and solving the full Navier-Stokes 
equations, the required mesh resolution is not practical. Thus, to maintain the stability 
and robustness of the numerical procedure, artificial dissipation is also included in viscous 
calculations. A semidiscrete form of (8) is then given by 

%- t (V-U) + Q(V)-D(U) = 0 (10) 

where V is the volume of the mesh cell being considered, Q is the discrete operator for the 
convective and physical diffusive terms, and D is the operator for the artificial dissipative 
terms. 

The artificial dissipation model considered in this paper is basically the one developed 
by Jameson, Schmidt, and Turkel [7]. This nonlinear model is a blending of second and 
fourth differences. The quantity D(U) in (10) is defined as 


D(U) = (D\ -D\ + Dl- D'JU 


( 11 ) 


where (£, 77) are arbitrary curvilinear coordinates, 




D\U = V S (A (+ • siJijjAtVeAjDy, 


( 12 ) 

(13) 


i,j are indices (for a cell center) associated with the £ and 77 directions, and A^, Vj are 
forward and backward difference operators in the £ direction. Following [6] and [8] the 
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variable scaling factor is defined as 


\+i,y - g + (*«)»+ wl 


(14) 


where _ 

(Af)t,; = ^«,;( r ) ’ 

M r ) - 1 + r ',i< 


(15) 


r is the ratio A,,/A £ , A £ and X n are the scaled spectral radii of the flux Jacobian matrices 
(associated with the £ and r? directions) for the Euler equations, and the exponent f is 
generally defined by | < f < 1. The spectral radii for the £ and rj directions are given by 


A{ = \uy n - vx n \ + cy/y* + x*, 


(16) 


A „ = |vx € - uy £ | + cy/xf+~yf, (17) 

and c is the speed of sound. The coefficients and use the pressure as a sensor for 
shocks and stagnation points, and they are defined as 



) max(i/,-_i ( y, Uij , i t y, i / »+ 2 1 j)’ 


(18) 
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Pi-1, i + 2p,-,y + Pi+l,;' ’ 


(19) 


(4) 

‘■+5. 


. = max 
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( 20 ) 


where typical values for the constants K^> and are j and y|g, respectively. 

For the normal direction ( 17 ), the dissipation contributions are defined in a similar way, 
except 


~~ *) * (Ar?)«,; 


( 21 ) 


Care must be taken, especially for internal flows, in defining the smoothing fluxes on the 
boundary. If the computation of the dissipative terms is carried out in each coordinate 
direction as the difference between first and third difference operators, conditions on those 
terms can be imposed at the boundary such that no errors in the conservation property 
are introduced, globally, by the numerical dissipation [9]. Moreover, such a scheme reduces 
the presence of undamped modes and guarantees a dissipative behavior for the dissipative 
operator D. 
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TIME-STEPPING SCHEME 

The system of differential equations of (10) are advanced in time towards the steady- 
state solution with a four stage Runge-Kutta scheme. This scheme is fourth order accurate 
in time for the linear problem and second order accurate for the nonlinear problem. It can 
be written for the time level n in the form 

u(°) = U n 

C/C 1 ) = IT® + a!AtR(UW) 

UW = Ef(°) + a 2 A tR{UW) 

( 22 ) 

17( s ) = Eft°> + ct z AtR(UW) 

UM = Ef(°) + a 4 AtR(UW) 

U n+1 = C7< 4 ) 

where at the ( q + l)st stage we have in the case of a single evaluation of the dissipation 

R{U^) = L{q(u^) - D{UW)), (23) 

and a good set for the coefficients a m (m = 1,4) is 

111 

<*i = -, a 2 = a 3 = -, a A = 1. (24) 

For economy the viscous terms are computed at the first stage and then frozen for the 
remaining stages. In multigrid computations, the driving scheme must be designed to 
rapidly damp out high frequency modes. Good high frequency damping properties are 
obtained with this scheme by evaluating the artificial dissipation terms on the first and 
second stages and freezing them for the remaining stages [10]. 

It is also important to note that the Runge-Kutta scheme has the desirable property 
that the steady-state solution is independent of the time step; therefore, the scheme is 
particularly amenable to convergence acceleration techniques. 

ACCELERATION TECHNIQUES 

Three methods are employed to accelerate convergence of the basic explicit time- 
stepping scheme. These techniques are as follows: 1) local time stepping; 2) residual 
smoothing; 3) multigrid. They are discussed in the subsequent subsections. 
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Local Time Stepping 

When the interest is only in the steady-state solution, a faster expulsion of disturbances 
can be achieved by locally varying the time step. In the present work the actual time step 
limit At is computed using 


At = cq 


/ At„ • A t d \ 
\ At c + A tj) 


(25) 


where A t c is the limit due to the convective terms, A t<j is the limit due to the diffusive 
terms, and Co is a constant usually taken to be the Courant-Friedrichs-Lewy (CFL) number. 
In particular, 



(26) 


At<j = 


V 2 


(27) 


K t fJ.[{x\ + y\) + {x* + y2)] ’ 
where A^, X n are defined in (16) and (17), respectively, and k t is a constant that has been 
set equal to 2.5 based on numerical experiments. 


Residual Smoothing 

The stability range of the basic time-stepping scheme can be extended using implicit 
smoothing of the residuals. This technique was first introduced by Lerat [11] for the Lax- 
Wendroff scheme, and later devised by Jameson [12] for the Runge-Kutta scheme. For 
two-dimensional flows the residual smoothing can be applied in the form 

(l — e^V^A^l — e^VqA^)R m = R m (28) 

where the residual R m is defined by 

Ri = ai^[Q(^)-L>(^], 

R m = a m ¥[Q(Ul m -')) - D(UW)}, (29) 

m = 2,4 

and computed in the Runge-Kutta stage m, and R m is the final residual at stage m after 
the sequence of smoothing in the £ and rj directions with the coefficients and e n . 

The use of constant coefficients in the implicit treatment has proven to be satisfactory 
(extending the CFL number by a factor of two to three) even for highly stretched meshes, 
provided additional dissipative support such as enthalpy damping [5] is introduced. How- 
ever, the use of enthalpy damping, which assumes constant total enthalpy throughout the 
flow field, precludes the solution of problems with heat transfer effects. The need for en- 
thalpy damping can be eliminated by using variable coefficients and e n that account 
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for the variation in cell aspect ratio. For the factorization of (28) effective expressions for 
the coefficients e(,e n can be derived by following the procedure of Martinelli in [6]. These 
expressions are written as 

(30)) 

where <j>(r) and <^(r -1 ) are the same quantities defined for the artificial dissipation, CFL is 
the local Courant number (usually taken to be 5), and the asterisk refers to the unsmoothed 
scheme. 

Multigrid Strategy 

Multigrid methods were first introduced for the solution of elliptic problems [13], and 
later Ni [14] and Jameson [10] applied them to the Euler equations. Actually, multigrid 
methods can be implemented to obtain improvement in the convergence rate in solving 
not only the Euler but also the Navier-Stokes equations [6,15]. The basic idea is to use 
coarser grids to speed up the propagation of the fine grid corrections. The Full Approxima- 
tion Storage (FAS) scheme in conjunction with the Runge-Kutta time stepping algorithm 
developed by Jameson [10] has proven to be an effective multigrid technique. 

For the multigrid process coarser meshes are obtained by doubling the mesh spacing. 
On the auxiliary meshes, the solution is initialized as 

„((>) _ W h U h 

U2h "no- 
where the subscript refers to the mesh spacing value and the sum is over the four fine grid 
cells that compose the 2 h grid cell. This rule conserves mass, momentum, and energy. On 
a coarse grid, a forcing function P is added to the governing discrete equations in order 
to impose the fine grid approximation. After the initialization of the coarse grid solution, 
this function is computed as follows: 

P2H = EMUh) - R2 h (U^). (32) 

Then, the time-stepping scheme on the ( q + l) at stage becomes 

pjf 11 = u2S - « s+ iA + P 2 ( ]?]. (33) 

We can also define a new value R* for the residual as 

R?h = P-2h{J^2h) + P2h, ( 34 ) 
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collect this value, restrict the solution to the next coarser grid, and repeat the process. 
The corrections computed on a coarse grid are transferred back to a finer grid with bilinear 
interpolation. 

Two types of simple fixed cycle multigrid strategies are primarily used. They are the 
V cycle and the W cycle. In the present work, results were obtained with both cycles. 
The structure of these cycles in the case of four grids is shown in Fig. (2). The W cycle is 
defined recursively since it becomes complex as the number of grids increases. At a given 
grid level subiterations (additional time steps) can be performed in both cycles. Moreover, 
the application of two Runge-Kutta (R-K) steps on the 2h grid and three R-K steps on all 
succeeding coarser grids is an effective strategy. 

In viscous flow calculations the viscous terms are computed also on the coarser grids, 
while the turbulent viscosity is evaluated only on the finest grid and then determined 
on each succeeding coarser grid by a simple averaging of surrounding finer grid values. 
The artificial dissipation model for the finest grid is replaced on the coarser grids with a 
simple constant coefficient second difference dissipation model. On each grid, the boundary 
conditions are computed in the same manner and updated at every Runge-Kutta stage. 
In order to provide an improved initial solution on the finest grid of the sequence of grids 
in the basic multigrid procedure, the Full Multigrid (FMG) method is employed. With 
the FMG method the solution is initialized on a coarser grid of the basic sequence of grids 
and iterated a prescribed number of cycles using the FAS scheme. The solution is then 
interpolated to the next finer grid. The process is repeated until the finest grid is reached. 
In the present work, three refinement levels are used. The first and second levels include 
two and three grids, respectively, and 50 cycles are performed on each. There are five grids 
in the final level. 

RESULTS AND DISCUSSION 
Bicircular Arc Cascade 

In [14] Ni presented results for inviscid flow over a 20 percent bicircular arc cascade 
with gap/chord ratio of 2. Then, Chima and Johnson [15] studied a 10 percent bicircular 
arc cascade for the case of laminar flow. The present numerical method has been used to 
study the first cascade geometry for laminar and turbulent flow conditions with separation. 
In Fig. 3(c) the 129 x 97 H-type grid used for the laminar calculation is presented. The 
spacing between the wall and the first cell center is 3 x 10 -4 chords and 10 -2 chords is the 
minimum spacing in the streamwise direction. For this case, the exit Reynolds number 
i?e 2 = 500, the exit isentropic Mach number (Af,) 2 = .5, and the angle of attack is zero 
degrees. The second difference artificial dissipation is set to zero. Figure 3(a) shows the 
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convergence history for the root mean square of the residual of the continuity equation 
scaled by the initial value. To indicate the monotone behavior of the multigrid scheme, 
300 iterations were performed on the finest grid. Using 5 grid levels and a V cycle with 
subiterations there are 7.2 decades decrease in the residual. In Fig. 3(b), the predicted 
isentropic Mach number distribution for the viscous solution is compared with that for 
the inviscid solution. The Euler solution was computed on a grid with the same number 
of cells but with a minimum spacing in each coordinate direction of 10 -2 chords. Mach 
number contours for this case are given in Fig. 3(d). The flow is separated at about the 
90 percent chord location. Due to the low Reynolds number, we found that this is one of 
the cases where it is important to compute the viscous terms on the coarser grids. 

Results for transonic turbulent flow (Mi = 0.675) are presented in Fig. 4. Figure 
4(c) displays the 129 x 97 mesh used in this test case. The minimum spacing at the wall 
boundary is 10 -4 chords, and in the x direction a spacing of about 1.5 x 10 -2 chords is used 
to obtain a good representation of the shock wave. For the inviscid mesh the minimum 
spacing at the wall is equal to the one in the x direction. Figure 4(a) shows the convergence 
history for a V cycle with subiterations. The variations in the computed isentropic Mach 
number are given in Fig. 4(b). In this case, the flow exhibits a little separation, involving 
about three cells, at the foot of the shock. Then, the flow reattaches and remains so until 
the 78 percent chord position. The Mach number contours for this flow (Fig. 4(d)) reveal 
a clean capturing of the shock. 

Gas Turbine Rotor Blade 

The present numerical method has been tested on several internal flow cases. In order 
to validate the method for turbomachinery applications, computations were performed for 
a Von Karman Institute (VKI) gas turbine rotor blade. This blade was tested extensively, 
as seen in [16,17]. At the inlet the flow angle with the respect to the turbine axes is 30 
degrees. The flow is turned about 96 degrees by the blade. For the solutions included 
herein the inlet Mach number is about 0.27 and the outlet isentropic Mach number (Mi ) 2 
varies from 0.81 to 1.20. Viscous results are compared with inviscid ones and with the 
available experimental data. 

The mesh used for the inviscid calculations is shown in Fig. 5. It is a 145 X 33 grid 
with minimum spacing in the x and y directions equal to 1.5 x 10 -2 chords. For the viscous 
computations (Fig. 6) the number of points in the y direction is increased to 65, and the 
spacing between the blade and the first cell center in the direction normal to the blade is 
set to 10 -4 chords. In both the meshes 32 points are located in the x direction before the 
blade, 65 points on the blade, and 48 points after the blade. A fine mesh is maintained 
after the blade in order to obtain good resolution of the shock system in the case with 
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supersonic outlet flow. 

In Fig. 7 the convergence history, Mach number distribution, and Mach number con- 
tours for the condition (Mi) 2 = .81 are displayed. The Navier-Stokes solution is in good 
agreement with the experimental data. This solution shows less tendency than the Euler 
one for overexpansion on the pressure side at the trailing edge. Some loss of accuracy in 
the trailing edge flow can be attributed to the mesh distortion of the H-type grid. Such 
distortion would also occur for a C-type grid. 

Results for transonic conditions ((M {) 2 = 1.0) are presented in Fig. 8. Again, there is 
good agreement between the viscous result and the experiment. Convergence histories for 
this case using the W cycle, V cycle with and without subiterations, and the basic scheme 
(CFL = 2.5) with only local time stepping are given in Fig. 9. Similar convergence behavior 
is obtained with both the W cycle and the V cycle with subiterations. Approximately 
the same amount of work (i.e., computer CPU time) is required for each cycle. The 
performance of the standard V cycle is not as good as that of the V cycle with subiterations. 
The convergence rate with each multigrid cycle is significantly faster than that of the basic 
scheme. 

In the case of supersonic outlet flow ((Af,-) 2 = 1.2) the numerical predictions are shown 
in Fig. 10. The agreement with the experimental measurements is good, and Mach number 
contours reveal a good resolution of the shock system. 

All the cases presented here suggest that engineering accuracy can be achieved in less 
than 100 multigrid cycles with the finest mesh. 

Supersonic Internal Flow 

Figures 11 and 12 show results for supersonic flow conditions. The geometry (see Fig. 
12) is the same as that investigated by Ni [14]; it is a duct with a 4 percent bump. No 
experimental data are available for this case, but comparisons between Euler solutions 
obtained from [14] and those from the present work indicate close agreement. Figures 
12(a) and (b) display respectively the viscous and inviscid Mach number contours for an 
inlet Mach number of 1.4. A 145 x 97 mesh was used for the viscous calculation, and a 
145 x 33 mesh was used for the inviscid one. For the viscous flow, — 1.5 x 10 ^ 

chords and Ay mt „ = 2 x 10 -4 chords, and for the inviscid flow, A z m ,- n = A y min = 1.5 x 10 -2 
chords. The two grids have a similar stretching in the y direction; this is why the shocks are 
smeared in almost the same way. To examine the shock capturing properties of the scheme, 
an inviscid computation on a weakly stretched 145 x 65 mesh was performed. As indicated 
in Fig. 12(c), the shock system is sharply captured with no tendency for smearing. Note 
that the same input coefficients were used in the dissipation model. Thus, care must be 
exercised in stretching the mesh when computing flows with strong oblique shocks. Next, 
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in Fig. 11(b) the Euler and Navier-Stokes solutions for the A/* distribution are presented. 
This disagreement is probably due to the interactions between shocks and boundary-layers. 
In fact, all the shocks cause separation at the locations where they impinge on the walls. 
The convergence history for the viscous case is given in Fig. 11(a). 

CONCLUSION 

The use of a Runge-Kutta scheme with acceleration techniques has proven to be fast, 
reliable, and accurate for the study of transonic internal flow. No problems lie in the 
extension of the numerical method to three-dimensional applications due to the versatility 
of the H-type grid. In addition, the very good convergence of the full multigrid process 
indicates the value of the method in turbomachinery component design. 
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Fig. 1: Auxiliary cell for computation of viscous terms. 
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Fig. 2: Schematic of V and W cycles with 4 grids. 
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Fig. 3: Bicircular arc cascade ((A/,-) 2 = .5 ;Re 2 = 500). 
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Fig. 5: 145 x 33 inviscid flow mesh for VKI LS59 gas turbine rotor blade. 
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Fig. 7: VKI LS59 turbine blade ((Aft) 2 = .81). 
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Fig. 8: VKI LS59 turbine blade ((iW,) 2 = 1). 
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Fig. 9: Convergence histories for VKI LS59 turbine blade [{Mi) 2 = 1). 
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Fig. 10: VKI LS59 turbine blade ([Mi ) 2 = 1.2). 
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Fig. 11: Supersonic internal flow (Mi = 1.4). 
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a) Viscous Mach number contours 



b) Inviscid Mach number contours 



c) Inviscid Mach number contours on weakly stretched mesh 


Fig. 12: Supersonic internal flow ( M\ = 1.4). 
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